Evolution to a smooth universe in an ekpyrotic contracting phase with w > 1 
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A period of slow contraction with equation of state w > 1, known as an ekpyrotic phase, has been 
shown to flatten and smooth the universe if it begins the phase with small perturbations. In this 
paper, we explore how robust and powerful the ekpyrotic smoothing mechanism is by beginning with 
highly inhomogeneous and anisotropic initial conditions and numerically solving for the subsequent 
evolution of the universe. Our studies, based on a universe with gravity plus a scalar field with a 
negative exponential potential, show that some regions become homogeneous and isotropic while 
others exhibit inhomogeneous and anisotropic behavior in which the scalar field behaves like a fluid 
with w — 1. We find that the ekpyrotic smoothing mechanism is robust in the sense that the ratio 
of the proper volume of the smooth to non-smooth region grows exponentially fast along time slices 
of constant mean curvature. 



I. INTRODUCTION 

For over two decades, the only known mechanism for 
homogenizing, isotropizing and flattening the universe 
was inflation, a period of accelerated expansion with an 
equation of state w (ratio of pressure to energy density) 
near -1. Its success in resolving the horizon and flatness 
problems is a principal reason why inflation became an 
essential part of the standard model of cosmology. In re- 
cent years, an alternative mechanism has been discovered 
in which smoothing and flattening occurs before the big 
bang as the universe undergoes a period of slow contrac- 
tion with w > 1. This alternative, known as the ekpyrotic 
mechanism [l|, Q , has been incorporated in alternatives 
to standard big bang inflationary cosmology including 
the "ekpyrotic" [l|, "new ekpyrotic" [|| and cyclic mod- 
els Q. We note that both the inflationary and ekpyrotic 
mechanisms can produce nearly scale-invariant spectra 
for density perturbations in addition to smoothing and 
flattening the universe. 

Until now, the ekpyrotic mechanism has only been 
shown to work in cases where the deviations from 
smoothness and flatness are small and perturbative when 
the ekpyrotic phase begins. The purpose of this paper is 
to show that the ekpyrotic mechanism is powerful and ro- 
bust enough to smooth the universe even when the initial 
perturbations are large and non-linear. 

Let us first review how the inflation and the ekpyrotic 
mechanism work in the perturbative regime where the 
cosmic evolution is well approximated by the Friedmann 
equation with an anisotropy term: 

3 V a a a 3(1+w ^ J a 2 a 6 w 



'Electronic address: garfinkl@oakland.edu 
^Electronic address: wlim@princeton.edu 
t Electronic address: fpretori@princeton.edu 
§ Electronic address: steinh@princeton.edu 



where H = a/ a is the Hubble parameter; a(t) is the 
Friedmann-Robertson- Walker scale factor normalized so 
that the value today, t , is a(t ) = 1; represents 
the present value of the energy density for component i, 
where m represents non-relativistic matter, r represents 
radiation and w represents an energy component with 
equation of state w, such as a scalar field and its poten- 
tial. The last two terms on the right-hand side represent 
spatial curvature and anisotropy. 

For an expanding universe, the term that dominates 
Eq. ([1]) after a long period of expansion is the one with 
the smallest power of a in the denominator. With only ra- 
diation and matter, the dominant term would be the spa- 
tial curvature, leading to a universe that is unacceptably 
open or closed by the present epoch. However, introduc- 
ing an energy component with w « — 1 totally changes 
the outcome because this component (p w ) then has the 
smallest exponent and dominates the Einstein equation, 
while the curvature and anisotropy (and other energy 
components) become negligible. This is the essence of 
how inflation works if the initial conditions are pertur- 
bative. For inflation, there are several "Cosmic No-hair" 
theorems in addition to numerical results supporting the 
claim that homogeneity, isotropy and flatness develop 
even when the initial conditions are non-linear and non- 
perturbative @, @. 

Now consider the analogous arguments for a contract- 
ing universe. With a(t) shrinking, the dominant term 
in Eq. (1) will be the one with the largest exponent of 
a in the denominator. For a universe with matter and 
radiation only, this would be the anisotropy term, which 
famously overtakes the evolution and drives the universe 
into chaotic mixmaster behavior. On the other hand, if 
there is an energy component with w > +1, then this en- 
ergy component dominates instead of anisotropy or spa- 
tial curvature, and chaotic mixmaster behavior never be- 
gins [IJ . For a scalar field 4> with potential energy density 
V((f>), the ratio of pressure to energy density is 
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which can be significantly greater than unity when V is 
less than zero and non-negligible and which approaches 
w = 1 if the scalar field kinetic energy dominates. The 
ekpyrotic phase in ekpyrotic and cyclic models includes 
an effective scalar field component of this type. 

For the cyclic model, the ekpyrotic phase is preceded 
by a period of dark energy domination and accelerated 
expansion which, if sustained long enough, would make 
the universe uniform and flat before the ekpyrotic con- 
traction phase begins. In this case, the perturbative ar- 
gument above should be reliable and sufficient to con- 
clude that the universe is smooth and flat as it approaches 
the big crunch. However, in the ekpyrotic or new ekpy- 
rotic models, generally, or in the cyclic model with a very 
short dark energy phase, the conditions at the beginning 
of the ekpyrotic phase are under less control. 

This paper investigates the robustness of the ekpyrotic 
smoothing and flattening mechanism when the initial 
conditions are non-linear and non- perturbative to deter- 
mine how the situation compares with the perturbative 
case and with inflation. For this study we are not con- 
cerned with any particular form for the initial conditions 
that may be motivated by some specific model; rather, 
we would like to understand how a generic, highly inho- 
mogeneous and anisotropic space-time evolves under the 
influence of the proposed smoothing mechanism, modeled 
here by a scalar field with a negative exponential poten- 
tial. Such negative exponential potentials arise naturally 
in supergravity and in string theory. Investigation of the 
proposed smoothing mechanism requires numerical solu- 
tion of the coupled Einstein-scalar system of equations. 
For simplicity, we restrict our studies to deviations from 
smoothness along a single spatial dimension. 

We use an orthonormal frame representation of the 
equations written in terms of Hubble-normalized, scale 
invariant variables 0] similar to that described in Q, 
though here coupling to a scalar field instead of a fluid, 
and using constant-mean-curvature (CMC) time slices. 
We discretize the equations using second-order accurate 
finite difference techniques, and solve them with a vari- 
ant of the Berger and Oliger [§] adaptive mesh refine- 
ment (AMR) algorithm for coupled elliptic-hyperbolic 
equations [1C| . We find smooth regions that are scalar 
field dominated in which the scalar field (kinetic plus po- 
tential energy density) component behaves like a fluid 
with w 3> 1, and also regions where the scalar field ki- 
netic energy dominates over the potential energy and the 
scalar field behaves like a fluid with w — 1. These lat- 
ter regions remain inhomogeneous and anisotropic, and 
throughout this paper we will refer to these parts of 
the universe as the "anisotropic regions" . Note how- 
ever that the anisotropic regions are neither anisotropy 
nor matter dominated because both the scalar field and 
the anisotropy of the metric play important roles in the 
dynamics. Futhermore, note that despite the fact that 
matter in the smooth regions behaves effectively like a 
fluid with io»l there is no issue of superluminal prop- 
agation as might arise from an actual fluid with such an 



equation of state. This is because we are always solv- 
ing the scalar wave equation with potential where dis- 
turbances always propagate within the light cone. In the 
anisotropic regions "spikes" also form, which are places 
where the fields change on very small spatial scales, and 
are similar to regions with this property that have been 
observerd in numerical simulations of singularities in vac- 
uum spacetimes [Hj . Despite the presence of the scalar 
field, the anisotropic regions exhibit dynamical behavior 
similar to chaotic mixmaster vacuum solutions, where 
there are a series of relatively quick transitions between 
longer epochs where the solution can be described by a 
w = 1 Bianchi type I spacetime. A difference here though 
is that there are only a finite number of transistions, so 
the mixmaster behavior terminates after several transi- 
tions. These dynamics are also known to occur in space- 
times where the matter is a fluid with w = 1 @, [HJ ■ AMR 
is necessary to resolve the spiky features that form both 
in the anisotropic regions and, in some instances, briefly 
in what will eventually become smooth scalar field dom- 
inated regions, and to resolve the almost domain wall- 
like transitions that develop between the smooth and 
anisotropic regions. 

The outline for the rest of the paper is as follows. In 
Sec. |II]we describe the equations, initial conditions, and 
numerical methods used to solve them. We present the 
results in Sec. IIIII The primary conclusion is that a 
scalar field with a potential inspired by cyclic models 
is a remarkably powerful and robust smoothing mecha- 
nism during a contracting phase of the universe, able to 
drive the spacetime to homogeneity and isostropy even 
starting with highly non-linear deviations from an FRW 
spacetime. Concluding remarks and a discussion of fu- 
ture work is given in Sec. IIVI 



II. THE EQUATIONS AND SOLUTION 
METHOD 

The method we use to evolve the Einstein-scalar equa- 
tions is the scale invariant tetrad method of Uggla et al. 
,7;]. We use this method with constant mean curvature 
slicing as is done in the vacuum simulations of [13l ] but 
with scalar field matter instead of vacuum. Thus our 
system can be thought of as the system of [l3[ but with 
extra variables and equations describing the matter, and 
with extra source terms for the influence of the matter 
on the metric evolution equations. More information on 
this type of method can be found in 0, [H[ . 

The spacetime is described in terms of a coordinate 
system (i, x l ) and a tetrad (eo, e Q ) where both the spatial 
coordinate index i and the spatial tetrad index a go from 
1 to 3. Choose eo to be hypersurface orthogonal with 
the relation between tetrad and coordinates of the form 
e = N~ 1 d t , and e a = e^di, where N is the lapse and 
the shift is chosen to be zero. Choose the spatial frame 
{e Q } to be Fermi propagated along the integral curves 
of eg. The commutators of the tetrad components are 
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decomposed as follows: 

[e ,e Q ] = u a e - (HS a fj + a a p )e p (3) 
[e a ,e^] = (2a[ a Sp] 7 + e af 3sn s ' 1 )e 1 , (4) 

where n af? is symmetric, and a al3 is symmetric and trace 
free. The scale invariant tetrad variables are defined by 
da = eo/ H and d a = e a /H while scale invariant versions 
of the other gravitational variables are given by 

{Ej,E a0 ,A a ,N a0 } = {e a l ,a a p,a a ,n af3 }/H. (5) 

Note that the relation between the scale invariant tetrad 
variables and the coordinate derivatives is 



d = N~ 1 d t 

TP if) 



(6) 
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where JV = NH is the scale invariant lapse. The matter 
model is a scalar field <j> with potential V of the form 



V(<P) = -V e 



-kip 



(8) 



where Vq and k are positive constants. The scale invari- 
ant matter variables are given by 



W = d cj> 
S a = d a (f> 
V = V/H 2 . 

The time coordinate t is chosen so that 

e -t = 3H. 



(9) 
(10) 
(11) 



(12) 



Here we have used the scale invariance of the physical 
system to make both t and H dimensionless quantities. 
Note that equation (|12|) means that the surfaces of con- 
stant time are constant mean curvature surfaces. Note 
also that the singularity is approached as t — > — oo. 



Due to equation (|12j) the scale invariant lapse satisfies an 
elliptic equation 

- d a d a N + 2A a d a N 

+ Af(3 + £ Q/3 £ Q/3 + W 2 - V) = 3. (13) 

The gravitational quantities E a \ A a ,N a @ and E a ^ sat- 
isfy the following hyperbolic evolution equations 



d t Ej = E a l ~M{E a l + i:jEp l ) 

d t A a = A, • |£ a 'd,.V d,,.V 

+ N (\d^J A a £ a % 

d t N af3 = N af3 - e^Zs^dyAf + Af(-N a/3 



(14) 
(15) 

(16) 



+ d <a d p> Af + A <a d p> N 



-3S 



a/3 



d <a A 



<aA3> 



- 2N <a ~ l N 0>1 + N~< 1 N <a p > 

+ e l5[a {d^N 0) s -2A^N fj) 5 ) + S <a S p> ] (17) 

Here parentheses around a pair of indices denote the sym- 
metric part, while angle brackets denote the symmetric 
trace-free part. 

The equations of motion for the matter variables are 
as follows: 



d t 4> = NW 
d t S a = S a • \\ 0,,X 



d a W - (S a + E Q %) 



(18) 
(19) 



d t W = W + S a d a Af 



dV 



+ N [d a S a -3W- 2A a S a - - . (20) 

In addition, the variables are subject to the vanishing of 
the following constraint quantities 



\ Xi 



e apx [d a E^ -AaE^-N^E^ (21) 

(Cj) 7 = d a N^ + t aPl d a Ap - 2A a N a ~ 1 (22) 

(C c ) a = d Ej - 3S Q % - e^N^s 1 

- WS a (23) 
C G = l + ld a A a -A a A a -\N afi N afi 

- \W 2 - \S a S a - \V (24) 
(Cs) a = S a -d a 4>. (25) 



The evolution equations can be freely modified by 
adding multiples of the constraints to them. In particu- 
lar, for numerical stability we add a multiple of (Cg) to 
the right hand side of the evolution equation for A a [1J| . 

The initial data must be chosen so that the constraint 
equations are satisfied, and then the evolution equations 
will ensure that they remain satisfied. We find solutions 
of the constraint equations essentially the same way as 
is done for more standard methods of numerical relativ- 
ity: by using the York method [l5j]. That is, we choose 
certain components of the variables and then solve the 
constraints for the rest. Our choice is by no means the 
most general possible one; but it is general enough that 
we expect any behavior that emerges from the evolution 
of these data to reflect the general behavior of singulari- 
ties for this type of matter. This expectation is bolstered 
by the experience of numerical simulations of vacuum sin- 
gularities. In particular, the greatest restriction on the 
generality of our initial data comes from the fact that we 
restrict our studies to deviations in homogeneity along 
a single spatial direction. Our spacetimes thus have two 
Killing fields. Nonetheless we choose our initial data to 
be the sort that in the vacuum two Killing field case[l6j 
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was sufficiently general that the behavior as the singular- 
ity was approached was the same as that of more general 
initial data in the case with no symmetries [13] • 

In the usual approach to numerical relativity the ini- 
tial data consists of the spatial metric and the extrinsic 
curvature. The York method then involves choosing a 
metric conformally related to the spatial metric and part 
of a tensor conformally related to the extrinsic curvature, 
and then solving for the conformal factor ip as well as the 
rest of the extrinsic curvature. Since we are using a tetrad 
approach, we must also have an initial spatial triad con- 
sistent with the initial spatial metric. For simplicity, we 
chose the initial conformal metric to be flat and (x, y, z) 
to be the usual cartesian coordinates for that metric, and 
we choose the spatial triad to lie along those spatial di- 
rections. Thus, the scale free spatial triad becomes 



It then follows from equation (0| that 
N af3 = 0. 



(26) 



(27) 
(28) 



The shear is essentially the trace-free part of the extrinsic 
curvature, and as in the usual approach in numerical rel- 
ativity, the constraint equations simplify for a particular 
rescaling of the trace-free part of the extrinsic curvature 
with the conformal factor. We therefore introduce the 
quantity Z a p defined by 



E a/3 = tp 



-6 7 
6 a/3- 



(29) 



Similar considerations apply to the matter variables, 
leading us to define the quantity Q given by 



w = V" 



(30) 



Here we will specify Q, <fi and a part of Zn~ and solve the 
constraint equations for the conformal factor ip and the 
rest of Zik- For convenience of the numerical simulations, 
we choose periodic boundary conditions < x < 2ir with 
and 2n identified where x is the single spatial coor- 
dinate that the metric and matter variables depend on. 
Also identifying the other spatial coordinates y and z 
means that our simulation is of a spacetime with spa- 
tial topology T 3 . Since the variables depend only on x 
and since x is periodically identified, specifying a vari- 
able means giving the coefficients of a Fourier expansion 
of that variable. 



From equation (|23|) and our ansatz for the scale invari- 
ant variables we obtain 



FZik = Qd k 



(31) 



In the vacuum case (i.e. for vanishing scalar field) 
this equation simply becomes the conditions that Z^ is 
divergence-free, which is in turn simply an algebraic con- 
dition on the Fourier coefficients of Z^. Note that since 



£ Q/ 3 must be trace- free, so must Z^.. A simple, but still 
fairly general divergence-free and trace-free Z^k is the fol- 
lowing: 



b2 £ 
'< ik = I £ a i cos x + bi 

g^cosx —bi 





O2 COS X 

- 62 — ai cos x , 



(32) 



where £, a\, 0,2, 61 and 62 are constants. We still keep 
this divergence-free part of Z^ but now add to it a piece 
that has a non-zero divergence. We simply specify the 
Fourier coefficients of <fi and Q via 



Q(x, t = 0) = -rjr cos(mia; + di) 
H 

(j)(x, t = 0) = /2 cos(m 2 x + d,2), 



(33) 
(34) 



where /1, mi, d\, /2, rri2 and di are constants. This turns 
equation (|31|) into an algebraic equation for the Fourier 
coefficients of this non-zero divergence piece of Z^ which 
we then solve. 

Now imposing equation (|24[) our ansatz yields 



z lk z lk )H 2 i>-\ 



(35) 



which is solved for the conformal factor tp using the nu- 
merical methods described below. 

The constraint equations (|2ip and (|22p are automati- 
cally satisfied by this ansatz. We then satisfy equation 
(I25p by using the given value of <fi to compute the initial 
value of S a . 



A. Numerical code 

We discretize the system of equations described in the 
previous section using second order accurate finite dif- 
ference methods, with Berger and Oliger [9fl style adap- 
tive mesh refinement (AMR) as provided by the PAMR 
toolkit [is) . On a single grid a two-time level, Crank- 
Nicholson(CN)-like discretization scheme is used. Stan- 
dard centered spatial derivative operators are employed, 
and in the hyperbolic evolution equations spatial deriva- 
tives (and undifferentiated functions) are averaged over 
the two time levels as usual within a CN scheme. Kreiss- 
Oliger dissipation [l!| is applied, and, although not nec- 
essary for the stability of unigrid evolutions, is important 
for reducing unphysical high-frequency solution compo- 
nents sometimes introduced at refinement boundaries. 
The elliptic equations are solved using an FAS (Full- 
Approximation- Storage) multigrid algorithm [20]. 

The elliptic slicing condition is incorporated into the 
Berger and Oliger time-stepping algorithm using the 
method described in [l(| • Such modifications are neces- 
sary to take advantage of time-subcycling, however here 
we find that we can evolve the system without time- 
subcycling yet keep the time step equal to that of the 
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coarsest level in the hierarchy. In other words, with a 
spatial refinement ratio of p sp a given level in the hi- 
erarchy (with level L\ being the coarsest) will have a 
CFL (Courant-Friedrichs-Lewy) factor Ai = Ati/Axi = 
Ati/Axi equal to p sv % ~ x times that of the base level Ai. 
In a typical simulation we use Ai = 0.2, p sp — 2, and 
have run cases where up to 25 levels of refinement were 
used, giving A25 ~3x 10 6 . Our code shows no signs of 
instability, and exhibits clear second order convergence. 
Though technically the evolution scheme is implicit due 
to the CN differencing, at each time step the code is 
able to converge to a solution within several iterations at 
most. We surmise that this rather atypical behavior for 
the solution of hyperbolic difference equations is due to 
the ultralocal nature of the spacetime in the approach to 
the singularity [ljl [I?) • This is reflected in the differen- 
tial equations by the spatial derivative terms becoming 
negligible, and hence are essentially reduced to a set of 
ordinary differential equations in time, one at each grid 
point in the domain. 




FIG. 1: t = const, snapshots of the normalized energy density 
in matter Q m (solid line), curvature fit (dot-dash line) and 
shear tt 3 (dashed line) for < x < 2ir at several times during 
the evolution of the initial data described in Sec lIIII Time 
runs from left to right along the top row and continues along 
the bottom row. The shaded slit (dotted outline) in the last 
panel (t — —150) indicates the range of x shown in the blow- 
up in Fig(2] 



III. RESULTS 



We have run simulations for a variety of initial condi- 
tions; here we show results from a single example that 
demonstrates the generic behavior : evolution from a 
highly inhomogeneous, anisotropic universe with signifi- 
cant curvature at the initial time to a universe containing 
distinct volumes of either smooth, homogeneous w 3> 1 
matter dominated regions, or w = 1 anisotropic regions. 
Whenever a w 3> 1 region forms it grows exponentially 
fast in proper volume relative to w = 1 regions. 

The particular initial conditions for this example are 



as 



«i 


= 0.70, 


«2 = 


0.10, £ = 


0.01, 


h 


= 1.80, 


6 2 = 


-0.15, 




h 


= 2.00, 


mi = 


=1, d!= 


-1.7, 


h 


= 0.15, 


m 2 = 


= 2, d 2 = 


-1.0, 



and 



V = 0.1, k = 10 



(36) 



for the scalar field potential parameters ([5]). The same 
initial data was evolved with several resolutions to con- 
firm second order convergence; the highest resolution has 
2049 points on the base level, and up to 12 additional lev- 
els of 2:1 refinement. 

It is enlightening to visualize the evolution via the be- 
havior of the matter (f2 m ), shear (£l s ) and curvature (Qk) 
contributions to the normalized energy density, defined 



-W 2 

6 vy 



1 na 



1 x^afl \ - 

~p a A a + A a A a 



(37) 
(38) 

(39) 



where Q, m +Q, s + Q,k = 1 by (|24|) . FigQ]shows these quan- 
tities plotted at select times during the evolution (which 
was stopped at t = —150). Note that all features in the 
solution are locally smooth — apparent step functions in 
some of the plots are simply due to the large size of the 
domain relative to the size of the feature. As an exam- 
ple, Fig[5]shows a zoom-in of the last panel of Fig[Tjabout 
one of the late-time spike structures that formed in the 
anisotropic regime. 

The effective equation of state parameter w is shown 
in FigO which takes the following form in Hubble nor- 
malized variables: 



jS a S a - V 



S a S a + V 



(40) 



By comparing FigQ] and Figj3] it is evident that at late 
times the region that has smoothed out and become 
matter dominated coincides with w 3> 1, whereas the 
anisotropic regime evolves to w — 1 (and exactly so to 
within numerical truncation error). 

The asymptotic behavior of the spacetime in the mat- 
ter dominated region appears to coincide with that of 
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X (blow-up of small interval) 

FIG. 2: Zooming in on one of the spike structures that formed 
in the anisotropic region at t = —150, as shown in Fig[T] 
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FIG. 3: The effective equation of state parameter w (|40[) for 
the simulation described in Sec lIIII for < x < 2ir at the 
same times as in Fig[T] At late times w — » fc 2 /3 — 1 in the 
matter dominated region, and w — ► 1 in the anisotropic region 
(in this simulation k = 10 for the potential (J8|). 



and the slicing condition for J\f (fT3]) becomes 

3 



3W ; 



3-V 



(42) 



where we have used f|41 1) to simplify the expression. We 
then assume that V remains finite and non-zero as the 
singularity is approached. This implies from (|8lllll2p 
that <f> takes the asymptotic form 



(f>{x,t) W (j) (x) 



2t 



(43) 



where k is the constant in the expression for the potential 
V = -V exp(~k(j)). Thus W © tends to 



W 



2 

W 



Combining these relations gives 



W 
V 



and from (|40p 



k, 

k 2 

3- — 
6 2 ' 

2 



fc 2 /3- 1. 



(44) 

(45) 
(46) 

(47) 
(48) 



We conjecture that at the singularity the above approx- 
imations become exact, namely 



lim W = 

t— oc 


k, 


(49) 


lim V = 


k 2 

3 - — 
6 2 ' 


(50) 


lim M = 

t — > — OQ 


2 

¥' 


(51) 


lim w — 


k 2 /3-l. 


(52) 



a spacetime with an isotropic singularity in the sense of 
Goode and Wainwright [2l| . We therefore conjecture that 
in an ekpyrotic phase with an exponential potential, an 
open set of initial conditions leads to an isotropic singu- 
larity. 

We can understand the behavior of the solution in the 
asymptotic matter dominated region by applying the fol- 
lowing approximation which we expect to hold to arbi- 
trary accuracy as the singularity is approached. To begin 
with, we neglect all spatial derivatives (such an approx- 
imation also holds in the anisotropic regions away from 
the isolated spikes). The constraint (|24[) then reduces to 



W 2 + 2V 



1«0, 



(41) 



Our simulations support this conjecture in that at late 
times (t = —150 in this example) these quantities have, to 
within numerical truncation error, reached their limiting 
values. Detailed analysis of the asymptotic dynamics can 
be carried out following the method of [12, 122, [23| . 



FigCH shows that as the singularity is approached the 
coordinate volumes of the it) = 1 vs. w 3> 1 regions of the 
universe are comparable. However, it turns out that the 
ratio of the proper volume of matter to anisotropic regions 
grows exponentially with time. Let S denote the proper 
spatial volume element associated with the spatial metric 
hij of t — const, slices, i.e., S = Vdet h. The fractional 
change of S with respect to time is 

1 



a In 5 



2 h l0 d t h l \ 



(53) 



7 



t=o 



1.5 
3N l 
0.5 

— 
1.5 - 
3N 1 - 
0.5 - 
- 



t = -7.5 



t = -2.5 



t = -10.0 



t = -5.0 



t = -150 



FIG. 4: 3 times the scale invariant lapse TV for the simulation 
described in Sec lIIII for < x < 2ty at the same output times 
as Figs Q]and 



which can be written as 

d t \nS = 3JV. 



(54) 
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Thus the scale invariant lapse Af is a direct measure of the 
rate at which the local proper volume element changes 
with time (and recall that t — > — oo as the singularity 
is approached) . FigJH shows 3Af from the simulation at 
several times. In the asymptotic regime where spatial 
gradients are negligible, N approaches a constant (|42|) . 
and thus (|54|) can be integrated to give 

S m oc e 6t/k \ w > 1 (55) 
S v ac e , w = 1 (56) 

where we have used (|42|) where w 3> 1, and note that 
V « when w = 1. Thus, at late times the ratio 1Z of 
the proper volume of matter to anisotropic regions of the 
universe grows as 

n= f&pte (xe -t { l-6/k>) (57) 

J S v dx 

Thus, as long as k > -\/6 (which is equivalent to w > 1), 
1Z — ► oo as t — > — oo. 

Fig[5] shows five state space orbits projected onto the 
(X + ,£_) plane, where 

S+ = i(Sn+S 2 2), E- = ^(Eu-E 22 ). (58) 

The orbits correspond to the evolution along the world- 
lines at x = 0,3.0,3.9,4.0,4.4. All five orbits begin in 
the upper left quadrant, away from the origin, indicat- 
ing anisotropic initial data. Towards the singularity, the 
first three orbits (solid lines) approach the origin of the 
plane, indicating isotropization. The fourth (dotted) and 
the fifth (dashed) do not isotropize. 



FIG. 5: The state space orbits for worldlines at x = 
0, 3.0, 3.9, 4.0, 4.4. Towards the singularity, the first three or- 
bits (solid lines) approach the origin of the plane, indicating 
isotropization. The fourth (dotted) and the fifth (dashed) do 
not isotropize. 



IV. CONCLUSIONS 

Our computations provide evidence that the ekpyrotic 
mechanism for smoothing and flattening the universe is 
robust and powerful, comparable qualitatively and quan- 
titatively to the inflationary mechanism incorporated in 
the conventional big bang model. This evidence is the be- 
havior as the singularity is approached of a class of space- 
times that, while not completely general, contain several 
degrees of freedom and begin far from FRW. Both the 
inflationary and the ekpyrotic mechanism require the ad- 
dition of an energy component that is commonly mocked 
up as a scalar field with potential energy. For inflation, 
the important feature is that, for some initial conditions, 
the scalar field can act like a fluid with w « — 1. It 
has been shown numerically that, beginning from highly 
non-linear and highly irregular initial conditions, regions 
dominated by the scalar field and with w « — 1 come 
to dominate the volume of the universe @. Similarly, we 
have found the regions in which the scalar field acts like 
a fluid with w > 1 come to dominate exponentially the 
volume of the universe during a contracting phase. 

This result addresses one of the key criticisms raised 
when the ekpyrotic model of the universe was first intro- 
duced; namely, it was suggested that the model required 
smooth initial conditions [24j . One of the motivations for 
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extending the ekpyrotic picture into a cyclic model was 
to include a period of dark energy domination before the 
ekpyrotic phase began in order to prepare smooth condi- 
tions Now, from the results here, it is clear that the 
dark energy epoch is not required for this purpose. Not 
only does this allow the possibility that the dark energy 
phase lasts only a few e-folds in the cyclic picture, as 
suggested in [2 51 ] , but it also opens the way for more gen- 
eral bouncing cosmologies that incorporate the ekpyrotic 
mechanism but do not cycle. 

With these results in hand, we are now prepared to 
tackle the bounce itself in the case that it is non-singular 
{a(t) shrinks to a non-zero value and then begins to in- 
crease). For the non-singular bounce, the equation of 
state must decrease from w > 1 to w < —1 for a fi- 
nite period during which anisotropy and inhomogeneity 



grows. Our goal is to determine if their growth can be 
kept at a level consistent with observations, establishing 
the viability of these bouncing cosmological models. 
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